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NOZZLE FLOW WITH VIBRATIONAL NONEQUILIBRIUM 


Introduction 

The research of this project concerns the modeling and numerical solution of the 
coupled system of compressible Navier-Stokes equations in cylindrical coordinates under 
conditions of equilibrium thermodynamics and nonequilibrium thermodynamics. 

The problem considered was the modeling a high temperature diatomic gas N 2 flowing 
through a converging-diverging high expansion nozzle. The problem was modeled in two 
ways. The first model uses a single temperature with variable specific heats as functions 
of this temperature. For the second model we assume that the various degrees of freedom 
all have a Boltzmann distribution and that there is a continuous redistribution of energy 
among the various degrees of freedom as the gas passes through the nozzle. Each degree of 
freedom is assumed to have its own temperature and consequently each system state can be 
characterized by these temperatures. This suggests the formulation of a second model with 
a vibrational degree of freedom along with a rotational-translation degree of freedom, each 
degree of freedom having its own temperature. Initially the vibrational degree of freedom 
is excited by heating the gas to a high temperature. As the high temperature gas passes 
through the nozzle throat there is a sudden drop in temperature along with a relaxation 
time for the vibrational degree of freedom to achieve equilibrium with the rotational- 
translation degree of freedom. That is, we assume that the temperature change upon 
passing through the throat is so great that the changes in the vibrational degree of freedom 
occur at a much slower pace and consequently lags behind the rotational-translational 
energy changes. This lag results in a finite relaxation time. In this context the term 
nonequilibrium is used to denote the fact that the energy content of the various degrees 
of freedom are characterized by two temperatures. We neglect any chemical reactions 
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which could also add nonequilibrium effects. We develop the energy equations for the 
nonequilibrium model from first principles. 

The basic equations describing the nozzle flow can be represented in various forms. 4 
This is done in order to check the derivations with other sources, references [2], [3]. The 
final form which is solved numerically are scaled equations in a dimensionless form. 

The two models to be compared are written in a weak conservative form using the 
symbols listed in the Appendix A. These derivations are from first principles and are con- 
sistent with other derivations as given in the references [18], [20]. The only difference in 
our derivation is that we use a more accurate representation for the relaxation time r as 
developed by Meador, references [7], [28]. The viscosity is obtained from standard formula- 
tions based upon the Sutherland potential, reference [13]. The thermal conductivities are 
developed by employing a Eucken approximation, reference [30]. We calculate the steady 
state solution based upon an assumed laminar boundary layer and do not define a Mach 
number because of the dispersive sound speed. Our objective is to determine the veloc- 
ity profiles and temperature profiles of both the translational-rotational and vibrational 
temperatures for the high expansion converging-diverging nozzle, which is defined in the 
Appendix B. 

The resulting equations which model the nozzle flow can be expressed in various forms 
as indicated in the following sections. In most forms the resulting equations are a coupled 
system of nonlinear partial differential equations subject to certain boundary conditions. 
To solve the resulting coupled system of nonlinear partial differential equations, several 
numerical techniques were investigated (i) The explicit MacCormack method, reference 
[14], (ii) The explicit-implicit MacCormack method, reference [14], (iii) The method of 
operator splitting, reference [14], (iv) Factorization schemes and (v) The Steger- Warming 
scheme, reference [25]. 

Single Temperature Equations 


For our first model we assume that there exists a single temperature T which charac- 
terizes the energy state of the system. The basic equations describing flow through a nozzle 
with cylindrical symmetry are given by (references [1], [3], [4], [12], [14], See also attached list 
of symbols in Appendix A.) 
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Continuity 

J( f e) + £(-flV r ) + £( r e v.) = o 

(1) 

Momentum 

~(r e Vr) + l(rT rr ) + §- 2 (rT rz )-T 6e =0 

4 

(2) 


at + Tt [ xTtz) + = 0 

(3) 

Energy 

J^(re t ) + j^(r[(« t + P)V r - V r r rr - V z T rz + 9r ]) 
+ ^(r[(e t + P)F* - V r Tj, r - V z t zz + q z \) = 0 

(4) 

where 

dT dT 

*— K SP ^— K Si 

(5) 


Trr = eV? +P-T rr 

(6) 


T rz — OYrVz T rz 

(7) 


Tee — P - Tee " 

(8) 


T Z z — + -P ~ T** 

(9) 

with the viscous stresses given by 



aVr 

T rr = + AV ■ V 

or 

(10) 


t„ = 2^ + XV-V 
dz 

(11) 


Tee = 2 v— + AV • V 
r 

(12) 


(dv z ^ av r \ 

r„ = r„ = ^ ar + 9g j 

(13) 


- ia, „ av. 

vv -; ft ( rv 'i + s ;- 

(14) 


The bulk viscosity is assumed to be zero, {Stoke ’s hypothesis), so that A = —2r)/3 and e* 
is the total energy per unit volume and is given by 


«t = * + f (V? + V z 2 ), (15) 

where e is the internal energy per unit mass determined from the relation 

de = C v dT. (16) 
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For a diatomic molecule the specific heats at constant volume and pressure are 


C v — C'urt and Cp — C%) d - P 


( 17 ) 


where R is the gas constant and C vrt = 5R/2 is the specific heat at constant volume due 
to rotational and translational degrees of freedom and 


Cim — R 



e HT 



(18) 


is the specific heat at constant volume due to the vibrational degree of freedom. The 
symbol (j> denotes the characteristic vibrational temperature which is unique for each gets 
species. For example, the characteristic vibrational temperatures for 0 2 is <t> = 2270° K 
and for iV 2 , <f> = 3395 0 K. For small temperatures C v « 5R/2 so that the vibrational 
degree of freedom only becomes excited when the temperature is on the order of <j>. 


Alternative Form for Energy Equation 

The energy equation has the Cartesian tensor form 

jj- + (e t V ( ) + ebiVi + (-PV { + mV,) ,i 


(19) 


where Q represents heat produced per unit volume, 6 t * represents body forces, and qi 
represents heat transfer per unit volume. The momentum equation is written 


dVi 

0 V ij V i ~ — P,j $ij + Tijj 


( 20 ) 


and when dotted with the velocity there results 


dV- 

qV~ + eViVjVij = QbiVi - ViPjSij + ViTiij 


( 21 ) 


or 


DV -* 

eV~ = e b-v-v-VP + v-v(T ii ). 


Consider the identity 


( 22 ) 


e~Me) = e 


dt 


(e t /e) + V • V(e t /o) 


de t e t dg -+ . . . 

= ar"7W + eF - v(et/e) 


( 23 ) 
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along with the continuity equation = — V(gV) and the identity 


V • (e t V) = V( e Ve t /e) = V(e t / e ) e ? + ^V(qV). 


These results show that 

e§- t (et/e) = ^ + ^V(eF) + W (e t V) - e jv( B v) 
= ^ + nVe t ) = ^ + (e t V i ), i 
Substituting equations (21) and (25) into equation (19) we obtain 


De D 


dQ 


D 


( 24 ) 


(25) 


eT^ + e^(V 2 /2) = '^--V-q+ S —(V i l2)+V.VP-V{PV) + {T ij V j ), i -V i r ijij . (26) 


Dt *Dt' dt 

Employing the identity 


PV -V = V(PV)-V -VP 
and defining the dissipation function 


^ — [*&) )t - ViUjj 


(27) 


the equation (26) simplifies to the form 

& + QV -Ve + PV-V + Vq = ?£- + $. 
at at 

When Q — 0 there results the alternative form of the energy equation 

e^ + f>F-Ve + PV-V + V-«-$ = 0, 

at 

where $ is the dissipation function given by 


(28) 


(29) 


* = V(r,-,^) - V • V(r iy ) = (r, y F y ), t - - V.n;,/ 


(30) 


and can be represented as 


$ = r/ [2 [D\ x + + ^ 33 ) + ( 2 D 12) 2 + (2JO13) 2 + (2P> 23 ) 2 ] + A0 2 (31) 


where 


Ay = -(V itj + V jti ) and 0 = At, 


(32) 
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is the dilatation. In the above equations the summation convention for repeated indices is 
to be understood. 

In cylindrical coordinates the dissipation function can be represented 

2 ' 


$ =77 

A 


<£>■}+(£ 


+ 


dV z 

dr 


+ 


In the limit as r — ► 0 we use L’Hospital’s rule that lim r _>o ^ 
The internal energy per unit mass is given by 


f T C v dT= f 
J T 0 Jt 0 


5R I R( i'T\* eHT 

— +R(4>/T) 


dT 


which integrates to 


1 Rd> de 

RT + ■ — dr Constant and — = C v . 


The various quantities in the energy equation (29) are given by 

V = V r e r + V z e z 

q — g r G r ~h Qz 

5 R<j> 

e = -RT + -r i T — 

2 e^' T - 1 

de _ de_dT __ dT 

dt ~ dT dt ~ ° v dt 

-* de de 

eV .Ve = e V r - + e V z - 


de_ r dT 
dr “ dr 


and 


de_ dT 
dz ” dz 


PV-V 
V •? = - 




av; r 

dr 

ld{rqr) , 3g 
r dr 


+ 


1 a , ^ , ar. 


The thermal conductivity K is a function of temperature T so that 

2 / om \ 2 


V • 9 = —if ( 


d 2 T ldT_ 
dr 2 + r dr + dz 2 y 


a 2 r \ 


aiiT 

dT 


dT 

dr 


) + (S 


(33) 


(34) 


(35) 


(36) 
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Vector Form For Equations of Motion 

The weak conservative form for the basic equations of motion can be written as 

w + L a J^l + ^ + L H > = 0 

dV r“ dr dz r“ 

where a = 0 is for two dimensional flow and a = 1 is for axisymmetric flow, and 


( 37 ) 


v' = coi(e, e Vr, ev z , e t ) 


(38) 


and 


G' = 


QVr 

gV r V r + P - T rr 
QV r V z - T rz 
L (e* + P)V r — V r T rr — V z T rz + q r J 


(39) 


with 


F' = 


qVm 


0 

ev r v z - r rz 

H’ = 

~P + Tee 

qV z V z + P -r zz 

0 

- (et + P)v z - V r Trz - V Z T ZZ + q z . 


0 


(40) 


The above equations are to be solved over the computational domain 0 < z < b and 
0 < r < f(z) where f(z) defines the shape of the nozzle. We introduce the change of 
variables 

x = z/b y = r/f(z), t = 




QobVp 

Po 


and write the system of equations in the form 


where 


dU dE dF 

~m + d^ + ^ + H - 0 ’ 


E=\f', F=jG'-^-F', H = f jF' + ±j{G‘ + H') 


for y 7^ 0. In the case y = 0 we use L’Hospital’s rule and find that 

__ . 1 ,dG' dH\ 

f * rdy + dy } 

We now introduce the scaled dimensionless variables 


u 1 


u 2 = 


QV r 


«3 


qVz et 

— — , u 4 = — 
0o Vo eto 


Qo QoV 0 * 

and then employ a numerical method to solve the resulting dimensionless system. 


(41) 


(42) 


(43) 


(44) 


(45) 
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Two Temperature Model 

For our second model we assume a vibrational degree of freedom together with a 
combined rotational-translational degree of freedom. Each degree of freedom is assumed 

A 

to follow a Boltzmann distribution and the energy content of each degree of freedom is 
characterized by temperatures T v and T respectively. As the gas passes through the nozzle 
there is a certain finite relaxation time r before the vibrational mode of excitation achieves 
equilibrium with the rotational-translational mode of excitation. Define the quantities: 


rii Population density of ith energy level 
Energy per molecule of the ith level 
h{ time rate of change of n t - due to V-T collisions 
qt Ordinary heat flux 


where q rr is the heat flux due to rotational energy 

q * Heat flux due to energy excitation of all energy levels 
* Total energy from all energy levels 

<f = 

i 

—♦ ^ ^ 

where U{ — V{ — V is the diffusion velocity of molecule in state % 

$ The Dissipation function. 


We construct a vibrational energy equation as follows. Let qe* — Yi n iU denote the total 
energy per unit volume from all excited states so that by integrating over the volume and 
surface of an arbitrary volume element we obtain 

J qe* dv = — j niCiVi • dS + n,c,* dv 


where dv is a volume element and dS is an area element of the control volume and n are 
rate equations to be determined. Using the Gauss divergence theorem and interchanging 
the order of summation and integration there results 


+ Y, V ( n * e ‘*(^ + V)) = Y hi€i 


^(fiO + VCee*?) = £ [n,fi - V(n,e;t?,)] . 


( 46 ) 
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Using the identities 


^(ee*) = ^§£1 + V ■ W) + V(ge*V) - ee*V - V — V ■ V( ee ‘) 

together with the continuity equation and 

D , De* >Dq De* - 

DV* } * Dt Dt * Dt K 

we write the vibrational energy equation as 

-Pe* V' * r? -* 

e pT = 2^ n * £ *- v ? 

t 

where the rate equations are from Meador, et. al. [28], are given by 

1 V ii r - e *~ e * 

t 

where the subscript e denotes equilibrium. In the case where * denotes the vibrational- 
vibrational energy mode (subscript v) we define 

e v{T) *- e v (T v ) DT V r _♦ n 

lim — — — - = C vv ——, lim q v = 0 

r— ►O T Dt r— >0 

so that the vibrational energy equation can be represented as 


De v DT V 

ur = eCvv ~Df 


e v (T) c w (Ty) 


-V-&. 


(47) 


The second energy equation is obtained by writing the energy equation (19) in the form 

ei-Aert + e„) + eV • V(« rt + e„) + i>V • 7 + V(f rt + <f„) - $ + eC„„X - eC m X = 0 (48) 
Ot 


where 


C,„,X 


e u (r v ) 


(49) 


and then employing the vibrational equation (47) to obtain the coupled energy equations 
- dert + e v • Ve rt + PV • V + V • q rt + $ + gC vv X =0 

( 50 ) 


dt 


de v 

lh 


+ qV * Ve v *t* V • q v — qC vv X — 0 
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which reduces to the form 

DT _ -1 
Dt gC rt 
DT V -1 


Dt qC t 


(v • q Tt + $ + PV • V + eC uv x) 
V • q v + X 


( 51 ). 


where 


= *(*) 


,4>/T 


(e*/T - 1)2 


<LT. 


(52) 


The integral in equation (52) is used to calculate 

1 


X = 


CvvT 


(e v (T) - e v (T v )). 


Integration produces the result 

T 2 1 — e - ^/ T * 


X = 


<j)T 1 — 


{“ p Ki - f )] - '} - 


(53) 


The quantity qC vv X = ^[e v (T) — e v (T v )) is thus a coupling term for energy between the 
vibrational and rotational-translational modes. The other terms in the coupled equations 
(51) are given by 


C v — C V rt Dvvi Dvrt — 5i2/2, 


0 e HT* 

c vv = R{<f>/T v ) 2 


e r t = -RT 


iv= L 


C,„, dT = 


R<j> 


To 

1 


e <(>/T* _ 1 


+ constant 


(54) 

(55) 

(56) 


X = -~(e v (T) - e v (T v )) 


X= fr (hiS) { 6XP [* (£ " ?)] - *} 


<frt = ~K rt VT 

<h = -x„vr u 

JT r t = 19rfk/4m 

Kv = rjCw 

ciff c r 3 / 2 


Sutherland’s formula 


' T + c 2 

C! = (1.488) * 2.16(10 _8 ), g c = 5J2/2, c 2 = 184.0 


(57) 

(58) 

(59) 

(60) 
(61) 

(62) 

(63) 
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Viscosity (kg/m-s) 



Temperature (K) 


Figure 1. Viscosity vs Temperature from Sutherland potential 


The Sutherland formula is illustrated in the figure 1 and depicts viscosity vs tempera- 
ture. Also in the above equations we have used the Eucken approximation, reference [30], 
that the coefficient of self diffusion D satisfies the relation pD = rj to obtain the specific 
heats K r t and K v . 
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The divergence of the heat flux terms are then given by 


V • Q r t — —K r t ^ 

V-g v = -K v (- 


dr 2 


dr 2 


iar 

d 2 T\ 

+ a*) 

dK rt 

r~dr 

dT 

ldT v 

, , d 2 T v 

\ dK v 

r dr 

^ dz 2 

) dT v 


(dT\ 
V Si- 


'S 2 (dT\ 

) + W 


(64) 

* 

(65) 


For the pressure we assume an equation of state for an ideal gas P = pRT. Following 
Meador et.al. [28] the relaxation time r for N 2 is taken as 


. 3.2188(10 -12 ) /T\ 1/2 , 

p{atm)r = I(T)Aw2T) \j) exp( " e/r) 


( 66 ) 


where <£,0, £ are characteristic temperatures given by <j> = 3395 K,0 — 3.2324(10 7 )K, 
£ = 95.9 K, and 

1/3 


t) 


1/2 


(l-e 2? ')exp[-(z + £+-C_)]cfo (67) 


where 




16 OTx 


L 27 <j>* 


(***)]■"• 


( 68 ) 


In the matrix form given by equation (42), we must replace the single temperature 
energy equation by energy equations associated with two temperatures. 

The energy equations can also be written in a conservative form, and are added to the 
continuity and momentum equations. The equations are then scaled and solved similar 
to the one temperature model. The conservative form of the two temperature energy 
equations are obtained by letting tt — E r t + E v and writing 


^Jr + + P)V r - VrTrr ~ V z T rz + g rt }+ 

d 1 

~q~ z K E rt + P ) V z ~ VrTrz “ V Z T ZZ + q zr t] + ~&C VV X = 0 

1 


where 


dE v 1 d , , d r 

-gf + ^d- r [rEvVr + + ai lE ’ v ’ + 9zv] 


E rt =- e RT+ |(V- r 2 + V, 2 ) 
E v = e iJ$/(e^ T ' - 1) 


e c vv x = 0 


(69) 

(70) 


(71) 

(72) 
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Both the single temperature and two temperature models have the weak conservative 
form given by the equations (42) . The single temperature model has column vectors of di- 
mension four, with while the two temperature model has column vectors of dimension five. 
For the two energy equation, nonequilibrium model, to facilitate the numerical method, 
the vector equations are written in the weak conservative form 


dUi 1 d(rGj) 
dt + r dr 


dFi 1 _ „ 

+ -r- 1 + -Hi = 0. 
dz r 


The continuity, momentum and energy equations then become 

dp 1 d(rpV r ) d(pV z ) _ 

dt r dr dz 


(73) 


(74) 


djpVr) 1 d{r(pV? + P- Trr)) , d(pV r V z - r„) 1 , p t Tee) _ 0 

dt r dr dz r ' ' 

d{pV z ) 1 d(r(pV r V z - r rz )) d(pV* + P-r„) 
dt r dr dz 

d{pe rt + f (V r 2 + Vj 8 )) 1 d(r((pe rt + f (V? + V?))V r - V M r„ + g rtr )) 

dt r dr 

. ^((pe„ -t- f (V^ -h ))V^ - ^ v n 

H ~ H pC vv X — 0 

dz 

^(P € v) , 1 d( r (P e vY 9ur)) , ^(P^v^z "h 9vz) v n 

1 1 — U. 

dt r dr dz 

The total equilibrium energy, E r t, is defined by 


(75) 

(76) 


(77) 

(78) 


E r t = pe r t + IjYr + Vz)i 


(79) 


and the total vibrational energy, E V) by 


E v pCy. 


(80) 


The resulting set of equations can be represented as a vector equation 


® .£!£) «: i 

dt r dr dz r 


(81) 


where U, F, G, and H are vectors. The vector U = col(/>, pV ri pV z1 E rt , E v ), is the set of 
conservative variables. 
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The equations axe now scaled by introducing the dimensionless variables defined by 



* _ A 

_ trt 

o 

~ Vq ’ 

Vr 

* _ e v 

_ Vug ’ 
L 

C U — y Q 2 » 

, v z 

~ Vo' 

Ql 

II 

# 

, T 

* n 
v = — • 

To’ 

Vo 


where 6 and L axe characteristic lengths, Vq is the characteristic velocity, po is the charac- 
teristic density, and To is a characteristic temperature. The characteristic viscosity, r) o, is 
the viscosity calculated at To. The Reynolds number, Re, is defined as the dimensionless 
combination of variables Re = foV ° L -. The nozzle domain becomes 0 < r* < ) and 

rjo o 

0 < z* < jj. The system of governing equations has the final scaled conservative form 




dt* 


dr * 


dz i 


r * 


where U*, G*, F *, and H* axe the column vectors 


U* = 


P' 

p*v; 

p*v; 

Kt 

L E* 


G* = 


P *V* 

P *v;v; + £p* - ££* 

n *v*v* — JLr* 

P v r v z HP R e 1 rz 


(K t + p*)v; - v; ir; - £v; ±t; z + g ; tr 


Etv;+g 


* 

vr 


F* = 


.* 

rz 


p*V* 

p*v;v; - ££«; 
p'v;v; + p* - Jjri 
(E; t + P')v; - v; i r ;, - v; + g ; tz 
eiv; + ?; 


i* 

fuar 


£ 2 p* 1 V 2 1 _ 


H* = 


0 


'P' Re'ee 


rL*p*CUX m 
rL* p*C* v X* 


(82) 


(83) 


(84) 


(85) 


( 86 ) 
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The scaled equation of state is 


P*=p*RT*( p|). 


(87) 


Finally, the nozzle is mapped to computational domain 0<®<l,0<y<l, which# 
is the unit square, by employing the change of variables, 


x = 


b/L’ V f(Lz*)/S’ 


( 88 ) 


The computational domain is then divided into 6 subregions for computational purposes. 
This computational domain is illustrated in the figure 2. The 6 subregions with nonuni- 
form grid spacing requires a weighted representation of the derivatives for computational 
purposes. The final system of equations has the computational form 


dU dE dF TT n 

- 4 - -}- — -}-* H — 0 

dt* dx dy 


where 


E — —F* F = — - — 
6 ’ f(bx) 


Lyf'(bx) 

G -~nb*r F ' 


H= w. r + 


(G* + JT), U = U\ 


f(bx) * ' yf(bx) 

When y approaches zero, we obtain by L’HospitaPs rule, the limit 


lim 


y— o 


yf(xb) 


(G* + JT) 


0 

L* 6*n 8*(i>V r ) 
pf^Re 5y 2 
__ L 2 6 djrrz) 

~S^ fRe dy 

L 2 6V, d(r rM ) L 3 K rt q 0 6 2 d 2 T ■ y 

6 2 fRe dy / 2 5y 2 ^ 

l? K v qp6 2 d 2 T v /~i y 

~~ £2 jl Qy2 (J\^ VV y\ 


(89) 


(90) 


(91) 
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Figure 2 . Computational Domain 
m 0 , mi, m2, m3 are number of divisions of x-region. 
no,ni,ri2 are number of divisions of y-region. 


Comparison with isentropic one-dimensional model 


For comparison purposes we also assume an isentropic process and calculate the results 
for a one-dimensional flow through the nozzle in the z direction where the area of the nozzle 
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is a function of z. For an isentropic process we have 

dh-d£- 

ds = = 0, 


with h = e + RT , e = C V (T) dT and dh = C p dT. Consequently, 

ds — C p ^ R p — 0 


or 


Let 


f F dP f T |7 /<A 2 e*/ T 

Jp 0 ^~Jto 2 + U/ (^ r - 


V - 


e ^/T _ 1 » 
and integrate by parts to obtain 


dF = 


1) 2 

e*! T <t> 

(e*/ T - 1)2 T* 


dr 
r ' 

dr 


_ 7 

IP„- 2 


iogP|?„=^iogr' T 


<t> 


1 ir + f t 1 dT 

r_x iTo +y To T 2 e 4,/T _ l dl - 


ITo re*/ 5. 

In the last integral, let z = e^/ r — l,dz = —t^l T -^dT so that 

f T 4 1 . r _ [‘ dz f Z — + f Z — 

Jto T 2 e4,/T ~ 1 Jzo 2 ( 2 + !) •/*<, 2 Jzo 2 + 1 

and consequently 


/* T 6 1 /l - e“^/ r ° \ 

J To Ti^HT^l dT = log ) • 


Therefore, we can calculate the pressure ratio as 


P /r\ 7/2 /l-e-*/ T « 
P 0 ~\T 0 ) \l-e-*/ T 

Since P = qRT we can write 


) (J&T ~ 


t/To 


e t/T 0 — 1 I ' 


Bo 


fT <t> \ 5/2 (l-e-*l T <‘\ ( 4>/T 4,/T 0 \ 

\<t>T 0 ) V 1 - e-*/ T ) “ P \e*/ T - 1 e^/ T ° — 1 / 


where <£/To is treated as a parameter. 

In one dimension the energy equation can be written 


dh + V z dV z — 0 or C p dT + V z dV z = 0. 


(92) 

(93) 

(94) 

(95) 


(96) 

(97) 

(98) 

(99) 
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Consequently, we may write 



V z dV z = - 



C p dT 


which integrates to 

- V 2 = 7R(To -T) + 2R4> ' ( 10 °) 


By dividing by a 2 = 'jRT, the local speed of sound, the one dimensional mach number 
can be represented 


M 2 = 


+ + L_ 

«y \ T J iT \e*/ T ° _ 1 e*/ T - 1 


( 101 ) 


with M = V-Ja. The mach number and one dimensional analysis is used to obtain an 
approximate solution to the more complicated two dimensional problem. Here 


C v 7 + 2(</>/T) 2 eV T (e*' T - l)" 2 
C„ “ 5 + 2(<I>/T) i e*/ T (e't’/ T - 1)“ 2 ‘ 


( 102 ) 


The one dimensional continuity equation is given by 


av z6 = a*v; 0 * 


(103) 


where the * quantities represent those values at the throat of the nozzle where M — 1 . 
That is, set M = 1 in equation (101), then solve the equations (101) (102) simultaneously 
for the value of (j>/T , treating <f>/To as a parameter. This calculated value of <j>/T gives 
T = T* when M = 1 and consequently we can calculate the values of P*, 0 *^*^V* at 
this critical value of the temperature. The equation (103) can then be expressed in the 
following form involving the above critical parameters 


A_ _ V z *g* _ - Rq*T* 

A* ~ V z q ~ v r \fiRT RqT^r 

y/lRT y/TRT* 


1 

L*t* p* i 

f*i*T P* 

(104) 

m\ 

/ lT pZjr~ M\ 

1 1 T* P 


1 

j'l* T <j> P* P 0 



m\ 

1 7 <t> T* Po P ' 
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Knowing the critical values P*, A% V* we can calculate the ratio T/4> as a function 

of A/ A* which is a function of z , with To/<t> as a parameter. These one dimensional values 
are then used as starting values (initial conditions) for the solution of the two dimensional 
non-isentropic nozzle problem. 

The figures 3,4 ,5 ,6 illustrate results from the above one-dimensional nozzle analysis. 


Boundary conditions 

Boundary conditions for both the model 1 and model 2 are similar. We state the 
boundary conditions for model 2. The steady solution of the governing equations is de- 
termined by the boundary conditions. In order to obtain a meaningful solution, boundary 
conditions must be applied that are not only applied correctly but are physically mean- 
ingful. The conditions are expressed in terms of the primitive variables [p,V r ,V X9 T,T v ,P\ 
and are converted to conservative variables when implemented numerically. 

At the inflow boundary, the density and the translational-rotational temperature is 
held fixed at some initial values po and To, respectively. This also means that the pressure is 
constant at the entrance. The vibrational temperature is assumed to be in equilibrium with 
the translational-rotational temperature and is also held constant at To. Gas is assumed to 
enter the nozzle parallel to the centerline, so the radial velocity is assumed to be zero. We 
set the flow to be subsonic at the inlet. An analysis of the flow characteristics of Euler’s 
equations, reference [27], suggests that one boundary condition must be left free to change 
with the solution of the interior flow. To this end, the axial velocity is extrapolated from 
interior data. In the computational domain, this represents the condition = 0. 

At the fax-field boundary, the flow is supersonic. All variables are extrapolated, or 


dp_ 

dx 


Since 


dV z _ dV r 
dx dx 


dP dT 

nR — 

dx P dx 


dT _ dT v _ 
dx dx 

(105) 

dp 

+ Er/ = o, 

dx 

(106) 


the pressure is also extrapolated. 

The symmetry of the nozzle is used to specify the conditions along the centerline. 
Assuming that the quantities above the center axis are mirrored by those below, we have 


dp _ dV r _ dV z _ dT _ dT v _ dP 
dr dr dr dr dr dr 


(107) 
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In the computational coordinates these conditions are simply 

dp _ dVr _ dV z dT _dT v _dP _ Q 

dy dy dy dy dy dy 

Furthermore, the radial velocity along the centerline is zero, 


(108) 


V r = 0 


( 109 ), 


since the radial velocity at the centerline must flow equally in the positive and negative 
directions. 

Lastly, boundary conditions on the nozzle wall are governed by the nozzle shape and 
the viscosity. Due to the viscosity, no-slip conditions are applied to the velocities at the 
wall boundaries, 

v r = v* = o. (no) 


The translational-rotational as well as the vibrational temperatures are extrapolated such 


that 


dT_-_ dT v _ Q 
dy dy 

The pressure gradient normal to the wall is assumed to be zero, 


dJP 

dn 


— 0 . 


(in) 


( 112 ) 


In computational coordinates this condition becomes 


dP dP [ 
dn dx 


(i +(m 


f(xb)f'(xb) 


dP 

dy 


0 . 


( 113 ) 


The density on the boundary is then calculated using the ideal gas law evaluated on the 
boundary. 

Initial conditions 

The choice of initial conditions is also important. Values that are too far from the 
steady solution can make the numerical solution unstable. And of course, the closer the 
initial conditions are to the steady solution, the less time is required for convergence. To 
minimize the convergence time, the initial conditions are the steady state solution values 
obtained from the one-dimensional model previously discussed. The initial vibrational 
temperature is assumed in equilibrium with the translational-rotational temperature before 
the throat. After the throat, the vibrational-temperature is frozen at the centerline, throat 
translational-rotational temperature. 
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In summary, the initial flow conditions were set along the centerline from a one- 
dimensional analysis of the subsonic to supersonic nozzle flow for the nozzle defined in 
Appendix B. These conditions were then extrapolated to fill the remaining nozzle grid 
points. A no slip velocity condition was applied at the nozzle walls and symmetry condi- 
tions were assumed along the nozzle centerline. The exit pressure was lowered to a point 
such that shockless supersonic flow was maintained in the diverging half of the nozzle. 
These boundary and initial conditions are consistent with other researchers, references 
[1], [2], [3), [9], [17], [20]. However, the choice of appropriately well posed boundary conditions 
for internal flows appears to still be an open question, reference [17]. 

Numerical Methods 

Various numerical methods can be used to solve the system given by the vector equa- 
tion (42) 


dU dE dF „ 

dt dx dy 

The following numerical techniques were applied to this equation. 


The explicit MacCormack method 

The explicit MacCormack method was developed in 1969 and is expressed by the 
algorithm: 

Predictor: 


U ? + 1 


XJV-. - 


At 

Ax 


(El 


i+l.J 


Eli) - - F Z ) - 


Corrector: 


l /. re + 1 = i 
>.-> 2 


where x t *,yy is the (*,j) node point. This scheme is second order accurate provided the 
various derivatives are differenced correctly. 
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Operator Splitting 


The weak conservative form for the equations (1) through (4) in terms of the compu- 
tational coordinates x, y are given by 


au 

dt 


+ 


dE dF TT n 

-r h -z 1- H — 0 

ox oy 


where 

E = F'/b, F = E’/b - F'yf'/f , H = H' + F’f'/f, 


with 


- co^rQ/roQo.rQVr/roQoVoyrgVz/roQoVoyret/roeto) 


TQVr 
To Co 

r(gv;?+P-T> r ) 
ro Co Vo 

y(gVrV»-T'r < ) 

Co Vo 

rUet + P)Vr-VrTrr-V,T rz +q r 


r -i 



r 0 Go 


r ° i 

r( e V r V,-r r ,) 


( — P+Too ) 

£ 

O 

CM 

H* ~ 

ro Co Vo 

"1 

+ 

1 

* 


0 

1 To Co Vo 

1 r((ct+P)V,-V r r r ,-V,r„+<j,) 


0 

rocto 




and 


E’ = 


We can then define the operators L x as the solution °f inf + ^ = 0 as 


Predictor: 

Corrector: 




i+1,3 


Ki) 


^ / - r , /j 


^ =2^y + U <J - ^ TO - 


Define the operator L y as the solution of ^ 4* = 0 as 


At 


Predictor: IT# =l£y - ~ Kj) 

Corrector: 


tt* * 




Define the operator L as the solution of ^ -j- = 0 as 

Predictor: =C// >y - AtH*j 

Corrector: U% =|(D? rf + AtF**) 
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where = F(Vtj),H% = H(U%), etc. 

The method of operator splitting requires that we time march according to the se- 
quence of operators 

U?f 2 =L x L y LLL y L x U? J : 


Factorization Schemes 

The system (42) is written as 


A U d_ 

At dx 





dH 

du 


AU = 0 


which can be written in the operator form 


a a 

I + At— A ■ +At-r-B ■ +AtC 
ox oy 


AU n+1 = —At AXJfj 


(114) 


(115) 


where A = £§, B = |C, C = are Jacobian matrices and 


AU ", = ( 


dE dF \ n 




Here we have used the notation 


(£-) 


A17 n+1 to mean — — ( AAU ) 
ox. 


For example, if A = A z is a forward difference operator we may write 

^.(AAU) » (AAU) = (A i+llj U i+hi - AijUij r)/Ax 


and if A = V x is a backward difference operator, we would write 

A(aAH) » ll (AAU) = (AijAUij - Ai-ijAUi-i^/Ax 

with similar results applicable to the forward and backward difference operators V y and 
A y in the y— direction. 

We desire to obtain the steady state solution to the system of equation (117). Toward 
this purpose, several factorization techniques have been investigated. Two such techniques 
are currently being used. One technique involves the predictor-corrector algorithm 
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Predictor: 


I + At^-A 
Ax 


+At-*-B-+AtC 

Ay 


At # 1 = -At AXJ”j 

vs'-iJTj+m 


n 

*»; 


where A x ,A y are forward difference operators in the x and y directions. 
Corrector: 


V- 

I + At-r-A 
Ax 


+At — ^£-+AfC 
Ay 

u n f x 


~ AU iJ 


=^K+c i +a c l ) 


where V x ,V y are backward difference operators in the x and y directions. This algorithm 
gives rise to upper and lower triangular systems of equations which can be solved at each 
time step. 

Another way to solve the system of equations (4) is to write it in the factored form 


7 . ,dA 
I + A*—— 

ox 


dB 

I + At— — I- AtC 
oy 


AU^f 1 - ~AtAU?<, 


Using central differences there results tridiagonal systems of equations to solve. Both of 
the above methods have been programmed and are working. However, these methods, as 
well as the previous methods, all suffer from the CFL (Courant-Fredrich-Lewy) time step 
restrictions which requires that small time steps be taken in order to achieve numerical 
stability. 

The explicit- implicit MacCormack method. 

The explicit-implicit predictor-corrector algorithm given by: 

Predictor with forward differencing 


/A 9 E* A y F* \ 

^ ***** = - At (-sf - + 


(116) 


with implicit system 


i 1 + At ^r + At ^r + Atc ) AU % +1 = (AC ^) 


Explicit 


(116a) 


with 


Z7*. n+1 = U£ + A C7*" +1 

ij »j * ij 


(1166) 
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du 

dy~ k 


h! + h 2 ( hi Ui ^ + h 2 ( UiJ 

1 / &2 / \ &1 , s \ 

ki + k 2 \ki ( UiJ+1 ~ Ui '^ + T 2 ^ Ui ' j ~ Ui ' j ~^ ) 


Figure 7. Irregular space grid points and first derivative approximations. 
Corrector with backward differencing 

(*£*+*£»♦«*») m 


with implicit system 


\ A V x -A A V y -B 
I + At--? + At—? 


+ A<cj AUV-t 1 = (A Dg.) 


Explicit 


tf" +1 = 5 [ds + (^')," +1 + a ^ +1 ] 


(1176) 
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The predictor corrector equations are applied to a rectangular t, j grid where 1 < i < 
ms and 1 < j < ns, with boundary terms given by i = 0,j = 0, i — msp = ms + l,j = 
nsp = ns+1. At certain %,j nodes of the grid we employ a variable grid spacing hi, j k 2 

a 

as illustrated in the figure 7 along with appropriate differencing relations. The equations 
(117), when expanded over the i,j grid, produces a system of equations to solve. In this 
system of equations there are boundary terms where along the boundary certain differences 
must be approximated, (reference [1]). 

Steger- Warming Vector Splitting 

The Steger- Warming Vector Splitting is used by Mr. J.G. Landry in his thesis to solve 
the nonequilibrium model 2. The method is described in his thesis. Upon completion of 
the thesis it will be added as an addendum to this report. 

Final form of governing equations 

In matrix form the basic equations are written as 

W IdW dF L 

dt' r“ dr dz r a 

where a = 0 is for two dimensional flow and a — 1 is for axisymmetric flow, and 

u = col( e, eV r , eV z , e t ) (us) 


and 




QVr 

f>V r V r + P-T rr 

0V rV z ~ T rz 

L (e* + P) V r - V r T rr - V z T rz + q r J 


with 


r qVz i 


0 

QV r V z + P-T rz 

H’ = 

—P + T ee 

qV z V z + P ~ t zz 

0 

- (e* + P)V Z — V r T rz — V z r zz + q z . 


0 


F* = 


We introduce the nondimens ional variables 
r 


r* = 


t* = 


V 


z * — ~ V* ~ Yl V — ~ n* = — 

_ L' r Vo’ z Vo’ P po 


t * rj 

L Vo rjo 


P* 


P + _ T * _ € 

T 0 ’ e = v$ 


then the above equations become 


dU* 1 d , dF * 1 „ 

dt* + r* a dr* ^ ^ + dz* + r* H ~ ’ 
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with 


U* = 


* 

a 


p'v; 1 

p’v; 

G* — 

p'v; 2 + P* - £r r * r 

p*v; 


p‘^7; - A*-;, 

< J 


( e t + nv; - -hv;^ - iv;t„ + g ; J 


and 


F* = 


— 1 
* 

<TS 

1 



0 

D *V*V* i-r* 

P Y r Y z Re rz 


1! 

* 

tel 

_ P* 4- J_ T * 

p’v;* + p* - ^4 

> 

0 

(«* + - JiV r *r;. - 7^*r* 2 + J 



0 


where 


«** = „v + ^-(v r * 2 + v; 2 ) 


„ .av,* x, / 1 a , dvt\ 
= 2ri 3 F + X + 

V v dr* + dz* j 


dF* 

r " = V dF + A * 


1 d . ± .. dv* 

-(r*v;) + ' 


r* dr' 


dx 


and 


0 

y* /I d dv*\ 

^ = 2,*^ + A*(--(r*V r *) + ^) 

* « *^v r * / 1 d . * Tr *. du!\ 

^ = V^ + V(- — (r*K r *) + ^J 


r) PoVoL * t7"2 \ \ > 

Fe = , et = etoe t) eto = Po Y o , A = j?oA 

Vo 


The shape of the nozzle is given by r = f(z) which becomes r* = f(bz*)/b in the * 
coordinates. Making the change of variable 


x—z 


V - 


f{bz*)/b' 


the final form of the equations are given by 


dU* d ,f> , by f . dF * b , TTt 6/' „ 

aF + &1 G ~ ~r F ) + ~di~ + Vf {H + G ) + T F = a 


In the limit as y — ► 0 we find that 

r 


lim 
v- o 


H* + G* 

yf 




dy 

0 


*tr* dV r 1_ & T rx 

z dy 


n*V 

P y * dy Re dy 

( e * p *) 9V r _ T rr 9V r ±-(V* ^ 

\^t^ r ) dy Re dv Re\ V i 


rx I a- -'*11 

Re dy Re\ T z dy 


+ Trz^fr) + 


By 


The above equations are to be solved over the computational domain 0 < z* < 1 and 
0 <r<f(z) where f(z) defines the shape of the nozzle. 


30 



SUMMARY MODEL 1 (EQUILIBRIUM THERMODYNAMICS) 


Continuity Equation 
Dp 


Dt 


+ p V • V = 0 


Equation of State 
P = pRT 


Momentum Equation 

DV 2 

p-=— = -V • P where P 1; = P6{j - + vy >t - - -&**>*,*) and rj is 

I/C 

Sutherland formula „ . 

ciSeT 3 / 2 , v . . 

*»= r + c ~ 

where for JV 2 we have c\ — 2.16(10~ 8 ) * 1.488, C 2 — 184.0 and g c = 32.174. 

Energy Equation 

= -[P : D + V • q- g rad ] 

where 

q = - AjfeVT 


_ 1 - . 

Dij — *** 

9rad =V (a R V(^ S bT 4 )) 




given by the 
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SUMMARY MODEL 2 (NONEQUILIBRIUM THERMODYNAMICS) 


Continuity Equation 

RE + ■ V = 0 

Dt 


Equation of State 
P = pRT 


Momentum Equation 

p— — = -V • P where P»y = PSij - T?(v»,y + vy,y - -6ikVk,k) and r? is 
Dt o 

Sutherland formula „ . 

e lSc T 3 / 2 llr , , 

1 s/m51 

where for AT 2 we have Ci = 2.16(10 -8 ) * 1.488, c 2 = 184.0 and g c = 32.174. 
Energy Equation 

D(C v T) 


and 


Dt 

DT V 

Dt 


= — [P : D + V • f r t — q ra d + pC vv X ] 

v.& 


rp2 

X=^~ 


<f> = 


+x 

1 — e~^/ r * 

<j>T 1 — e~^/ T 

hu 




=fc exp(^-) ( exp ^^ _ 1 

Qv — — A v vr v 
Ay —TjC vv /m 
qn = - AjtVT 


-2 


Afc =19rjk/4m 
1 

mw p(atm) 


exp{A(T~ 1/3 - 0.015 /z 1/4 ) - 18.42} 


CO 

to 

1 — exp(— ^/5000) 

V5oooy 

. l-exp(-(£/r) . 


A =220, <j) = 3395, ^ = 14. 


Ay ^Ky + «y ,<) 


Qrad 


=V (^V^assT 4 )) 




given by the 
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Preliminary Considerations 

The Knudsen number, reference [31], is defined K n — A/d, and is the ratio of the mean 
free path A to some typical dimension in the flow field. The Navier-Stokes equations, for 
the simulation of gas flows, is valid when the Knudsen number is very small in comparison 
with unity. As the Knudsen number increases the intermolecular collisions of the gas 
particles becomes less. A knudsen number of 0.1 has been used, reference [29], to describe 
the boundary between continuum and transition regime flows. This assumption is based 
upon the selection of an appropriate scale length d used to determine the Knudsen number. 

For internal flows, gas properties which depend upon molecular collisions such as 
viscosity, heat conduction, diffusion, heat capacity, etc., will be greatly altered when the 
Knudsen number is large. This is because molecules will collide more readily with the 
wall boundaries than with each other. As the density of the gas decreases, the collision 
rate of the molecules diminishes and eventually the continuum theory breaks down. Under 
such conditions one can resort to direct simulation Monte Carlo techniques associated with 
Boltzmann equation, as opposed to the Navier-Stokes equations, which treats the rarefied 
gas as a set of discrete molecules. 

The maximum mean free path for N 2 is given approximately by 

A max (cm) ^ o 

reference [26], where a 2 = 14.9 (10) ~ 16 cm 2 and n is the gas density in moleculues/cm 3 . 
The figure 8 illustrates the maximum mean free path vs density of Nitrogen. Using the 
local nozzle minimum radius as representative of d, and using the density from a one 
dimensional nozzle analysis, the figures 9,10,11 illustrate approximate Knudsen numbers 
vs distance along the nozzle. These figures illustrate that for high expansion nozzles, large 
pressures are required to insure the Navier-Stokes modeling remains valid. 
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Conclusions 

The reference [l] points out several difficulties associated with the numerical solution 
of the compressible Navier-Stokes equations. These difficulties are: 

(i) The Courant-Friedrich-Lewy (CFL) condition for numerical stability limits the time 
step that the solution can be advanced. 

(ii) The treatment of steep flow gradients and mesh size needed to handle these gradients. 

(iii) The treatment of turbulence and associated mesh size. 

(iv) Approximate factorization error in the numerical solution. 

(v) Linearization error. 

The figure 2 illustrates the mesh selected for the numerical solution in order to address 
the concerns of (ii) and (iii) above. The size of the mesh was selected to address the 
concerns of items (iv) and (v) above. The step size limit continued to be a problem in our 
numerical techniques for a solution. 

The diffusion fluxes, described by the vectors E, F and H of equations (42) (43), are 
found from the Navier-Stokes diffusion of momentum and Fourier’s heat conduction law. 
The viscous terms and heat conductivity, in these terms, are all temperature dependent. 
The heat capacity of the gas is taken from reference [19]. We also employ the Stokes 
hypothesis that the second coefficient of viscosity is given by A = —2*7/3. The exit pressure 
is such that a shockless supersonic flow is sustained in the diverging portion of the nozzle. 

The reference [17] points out that there are no proper initial and boundary conditions 
that will insure existence and uniqueness of the solution. For this problem we tried to 
select boundary conditions that were physically realizable, we use adiabatic boundary 
conditions for the temperature T and required that the normal derivative of the pressure 
be zero at the walls. The exit values on the boundary were all extrapolated. The inlet 
boundary conditions were such that all values were fixed except for one variable conditions 
which was extrapolated (reference [27]). We investigated two different input conditions. 
One the extrapolation on density and the other investigated the extrapolation of velocity. 
The density extrapolation seemed to keep the mass flow, (continuity equation), in balance. 

For model 1 we employed the explicit-implicit MacCormack method. For the model 
2, the Steger- Warming technique was employed. The numerical results for the model 1 
are summarized in the figures 12 through 19. The figures 12,13 illustrate the resulting 
temperature contours throughout the nozzle. Note the sudden drop in temperature as 
the gas passes through the throat. The figures 14,15 illustrate the contour plots of the 
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resulting radial velocity V r , while figures 16,17 illustrate the contour plots of the axial 
velocity V z . The figures 18,19 illustrate the contour plots of the logarithm of density. All 
nozzle dimensions are given in the appendix B. 
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r Contour Plot of Temperature 
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Figure 12. Contour plot of temperature. 
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Figure 13. Contour plot of temperature. 
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Figure 18 . Contour plot of logarithm of density. 
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Figure 19. Contour plot of logarithm of density. 
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The model 2 assumes nonequilibrium thermodynamics. The vibrational potential of 
the molecules are represented by a harmonic oscillator model, references [21], [28]. The 
nonequilibrium thermodynamic model is formulated with two temperatures and differs 
from the models developed in the references [2] and [3]. The main differences are (i) ther- 
mochemical vs thermodynamic nonequilibrium, (ii) boundary conditions, (isothermal vs 
adiabatic.) (iii) nozzle shape, scaling and mesh size, (iv) The representation of relaxation 
time over the solution domain. 

The nonequilibrium model is developed from first principles. It is a two tempera- 
ture model used to study thermodynamic nonequilibrium where vibrational excitation and 
rotational-translational excitation are the dominant nonequilibrium phenomena. Such sit- 
uations arise in the study of re-entry flow and certain converging-diverging nozzle flows, 
reference [2]. Almost everything in model 2 is the same as in model 1 with the exception 
of energy and temperature representation. The vibrational relaxation follows the work of 
Meador, et. al., reference [7], for Nitrogen gas through a nozzle. 

The numerical results for model 2 are summarized in the figures 20 through 29. The 
figures 20,21 are contour plots of temperature. The figures 22,23 are contour plots of 
vibrational temperature. The figures 24,25 are contour plots of radial velocity V r . The 
figures 26,27 are contour plots of the axial velocity V z . The figures 28,29 are contour plots 
of logarithm of density. 
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Figure 20. Contour plot of rotational-translational temperature. 
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Figure 21. Contour plot of rotational-translational temperature. 
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Figure 22. Contour plot of vibrational temperature. 
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Figure 23. Contour plot of vibrational temperature. 























Results indicate that the vibrational temperature has essentially a frozen value, along 
the centerline, upon passage through the throat area into the diverging portion of the 
nozzle, the same as reported in the references [2], [3]. As the flow expands into the diverging 

s 

nozzle region, the density of the gas decreases, and the Knudsen number becomes large. 
Thus, care must be exercised in the nozzle design (Shape of diverging portion of the nozzle) 
and choice of initial parameters (Inlet pressure) in order that the model remain valid. 

In the boundary layer the vibrational temperature is not frozen. These vibrational 
temperature effects imply that there is no meaning to the term ”Mach Number” in this 
situation as the speed of sound is dependent upon the vibrational properties of the gas. 
In both models the boundary layer thickness becomes very thin for high Reynolds number 
flow. A high initial pressure insures both a low Knudsen number and a thin boundary 
layer. Also, an increased pressure produces a more uniform flow in the nozzle. These 
results are consistent with those reported in reference [5]. 

A more detailed analysis of model 2 will be given in Mr. Landry’s thesis, which will 
be an addendum to this report and submitted at a later date. 

In conclusion, we have developed two models for the study of high pressure gas flow 
through nozzles. We have investigated the basic physics associated with the modeling of 
such a problem under a variety of circumstances. Improved models for gas vibrational 
relaxation were employed during the analysis. The numerical analysis produced both 
quantitative and qualitative results associated with both equilibrium and nonequilibrium 
thermodynamic flow physics through a converging-diverging nozzle. It turns out that these 
problems are extremely difficult to solve numerically and care should be taken upon the 
selection of the numerical method employed. They are challenging problems for current 
and existing numerical techniques. The development of numerical techniques which will 
handle these types of problems quickly and efficiently are areas for additional research. 
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APPENDIX A 
LIST OF SYMBOLS 


a 

A 

b 

0 V1 C vr t , Vyy 

A,; 

D d - „ 
— = — + Y • V 

d* a* 

€ rt j € u 


h 

h 

k 

K,K rU K v 
m = VF/iV 0 
M 
Na 

Qi Qrty <7v 


P 

P 

r 


t 


T,T V 


Speed of Sound [m/s] 

Cross sectional area [m 2 ] 

Body force per unit mass [Newton/ Kg] 

Specific heat at constant volume [Joule/ Kg K] 
Rate of deformation tensor [s _1 ] 

Material or substantial derivative 
Energy per unit mass [Joule/ Kg] 

Energy per unit volume [Joule/m 3 ] 

Plank’s constant [Joules] 

Enthalpy [J oul e / m 3 s] 

Boltzmann’s constant [Joule/ K] 

Thermal conductivities [W/mK] 

Molecular mass [Kg] 

Mach number 

Avogadro’s number [mol -1 ] 

Heat input per unit volume [Joule/m 3 s] 
Pressure [Newton /m 2 ] 

Gas Constant [Joule/ Kg K] 

Radial distance [m] 

Entropy per unit volume [Joule/m 3 s K] 

Time [s] 

Temperatures [K] 


A 
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V Velocity [m/s] 

V r ,V z Velocity components [m/s] 

W Molecular weight of N2 [Kg] 

X Coupling term [K/s] 
x, y Computational coordinates 
z axial distance [m] 
r 1 Viscosity coefficient [Kg /ms] 

X Second coefficient of viscosity [Kg/ 
g Density [Kg/m 2 ] 
t Relaxation time [s] 

Tij Stress tensor [Newton /m 2 ] 

(f), 0 , £ Characteristic temperatures [K] 
v Frequency [s _1 ] 

__ 

C v 


Ratio of specific heats 



APPENDIX B 
NOZZLE DESIGN 

The nozzle shape is illustrated in the figure B1 and can be described by a straight 

A 

line 1 , circle, cubic spline and straight line 2. In this figure a 0 = 0.5 mm is the throat 
radius, 6 = 350mm is the length of the nozzle, ro = 5.0mm is the entrance radius, and 
R = 2.0mm is the radius of the circle. The parameters di,d2,d3,d4 a-re selected such that 
the line 1, circle, cubic spline and line 2 are connected in a continuous fashion. We require 
that at z = ai that the slope of the line 1 and slope of the circle are equal. At z — d2 the 
slope of the circle is zero. The parameter 03 is selected such that the cubic spline converts 
the circle to the line 2. 

We write the equation of the circle as 

(r - (a 0 + R )) 2 + (z - a2) 2 = R 2 (61) 

and the equation of the line 1 as 

r — 5 = ml * z (52) 

where mi = — tan(60°) is the slope of line 1. The slope at any point on the circle can be 
determined from the derivative of equation (bl) 

(r — (ao + R))— + [z — 02) = 0 (63) 

At z — d2 we require that the slope of the circle is zero, = 0 and at z = ai we require 
that the point z = «i,r = ri lie on both the circle and line 1 so that 

ri - (a 0 + R) = —y/R 2 - (di - d 0 ) 2 (64) 

and 

\JR 2 — (di — d2) 2 * mi = — (di — d2). (65) 

Solving the equation (b4) for di — d2 we find that 

di — d2 = —VzR/2. (66) 


From the equation (b4) we find 



( 68 ) 


The circle intersects the line 1 at the point (ai,ri) which gives the equation of line 1 

r - (a 0 + R - y/ R 2 - 3) — mi * (z - ai). 

Note that at z — 0, where r = ro = 5 we find that 

5 — (ao + R ~ y/r 2 — Z) = y/Zai 


or a i 


5 — (ao + R — yj R 2 — 3) 


V5 


For 0 < z < a\ the nozzle shape is described 


(69) 


r = f[z) — mi * (z — ai) + ao -f- R — yj R 2 — 3 


(610) 


or . 

r = f(z) = mi * (z — 02) + ao + R ~ 3 — y/R 2 — 3 (611) 

where mi = -tan(60°) = -\/3. For the region d\<z<d 2 the nozzle shape is described 

r = f(z) — do + R ~ y/R 2 — (z — d 2 ) 2 (612) 

For the region 02 < z < 03 we must find the parameter 03 such that the cubic spline 
converts the circle to the line 2. The cubic spline is represented els 


r = £3(2) = A(z - a 2 ) 3 + £(z - a 2 ) 2 + C(z - a 2 ) + D 


(613) 


and is subject to the end conditions 


dt 

£3(02) = ^U=o 3 =0 


£3(^2) ^0 £3(03) =a4 (unknown) 


£"(a 2 ) = 


R 


£3(03) =m 2 = tan0 (slope of line 2) 
£3* (^3) = o 


Solving this system of equations we find 

-1 


A = — —r B — — C — 0 D — a 0 

12R 2 m 2 2R 


(614) 


and consequently 


dz= d 2 -\- 2 Rm 2i a 4 = £3(03) = ao + 4Rm\/Z 


(615) 


63 



Then the nozzle shape for the region a 2 < z < a 3 is given by 


-(z-a 2 ) 3 (z - a 2 ) 2 
r - f(z) = V + <*o 


12# 2 m 2 


2i2 


(616) 


and for the region a 3 < 2 < 6 we have 


r = f(z) = a 4 + m 2 * (z - a 3 ). 


(617) 


The equations (bl0),(bl2),(bl6) and (bl7) must be converted from millimeters to meters 
for use in our analysis. 

In summary, for r in meters, z in millimeters and using the conversion factor cv = 
10~ 3 , the equations describing the shape of the nozzle can be written 


0 < z < ai, 
a\ < z < a 2i 


a 2 < z < a 3 , 
a 3 < z < 6, 


r = f(z) = cv * (m 1 * (z — a 2 ) — a 0 + i? — 3 — \/R 2 — 3) 

r = /(*) = cv * (a 0 + JR — y/R 2 — (z — a 2 ) 2 ) 

, v , — (2 — a 2 ) 3 (2 — a 2 ) 2 \ 

r=/w=ro ‘W + ir +ao) 

r = 7(2) = cv * (04 + m 2 * (2 — a 3 )) 


where 


mi — — tan 60° 
m 2 = tan 0 
6 =4° 

a,Q =0.5 mm 
R =2.0 mm 
ro =5.0 mm 


01 


5 — clq — R + y/ J2 2 — 3 


o 2 =di -j- \/3.R/2 


a 3 =a 2 + 2i2m 2 
ARm/k 

a 4 =a o H — 
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